{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Classification example for Landsat 8 imagery\n",
    "\n",
    "**Author:** René Kopeinig<br>\n",
    "**Description:** Classification example for Landsat 8 imagery based on the scientfic work \"MAD-MEX: Automatic Wall-to-Wall Land Cover Monitoring for the Mexican REDD-MRV Program Using All Landsat Data\" by S.Gebhardt et. al 2014. Please find the link to the paper here: https://www.mdpi.com/2072-4292/6/5/3923\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import Image\n",
    "import ee\n",
    "ee.Initialize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get and visualize the Landsat 8 input data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "area_of_interest = ee.Geometry.Rectangle([-98.75, 19.15, -98.15,18.75])\n",
    "mexico_landcover_2010_landsat = ee.Image(\"users/renekope/MEX_LC_2010_Landsat_v43\").clip(area_of_interest)\n",
    "landsat8_collection = ee.ImageCollection('LANDSAT/LC8_L1T_TOA').filterDate('2016-01-01', '2018-04-19').min()\n",
    "landsat8_collection = landsat8_collection.slice(0,9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"https://earthengine.googleapis.com/api/thumb?thumbid=d1d1930789549e3839ca2eb086bf521e&token=0bc2ded045d0168feff6a6961a18440b\"/>"
      ],
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vis = {\n",
    "    'bands': ['B6', 'B5', 'B2'],\n",
    "    'min': 0,\n",
    "    'max': 0.5,\n",
    "    'gamma': [0.95, 1.1, 1],\n",
    "    'region':area_of_interest.getInfo()['coordinates']} \n",
    "image = landsat8_collection.clip(area_of_interest)\n",
    "theGeom = area_of_interest.getInfo()['coordinates']\n",
    "thumbnail = image.getThumbUrl(vis)\n",
    "Image(url=thumbnail)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Functions to derive vegetation indices and other raster operations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def NDVI(image):\n",
    "    return image.normalizedDifference(['B5', 'B4'])\n",
    "\n",
    "def SAM(image):\n",
    "    band1 = image.select(\"B1\")\n",
    "    bandn = image.select(\"B2\",\"B3\",\"B4\",\"B5\",\"B6\",\"B7\",\"B8\",\"B9\");\n",
    "    maxObjSize = 256;\n",
    "    b = band1.divide(bandn);\n",
    "    spectralAngleMap = b.atan();\n",
    "    spectralAngleMap_sin = spectralAngleMap.sin();\n",
    "    spectralAngleMap_cos = spectralAngleMap.cos();\n",
    "    sum_cos = spectralAngleMap_cos.reduce(ee.call(\"Reducer.sum\"));\n",
    "    sum_sin = spectralAngleMap_sin.reduce(ee.call(\"Reducer.sum\"));\n",
    "    return ee.Image.cat(sum_sin, sum_cos, spectralAngleMap_sin, spectralAngleMap_cos);\n",
    "\n",
    "#Enhanced Vegetation Index\n",
    "def EVI(image):\n",
    "    # L(Canopy background)\n",
    "    # C1,C2(Coefficients of aerosol resistance term)\n",
    "    # GainFactor(Gain or scaling factor)\n",
    "    gain_factor = ee.Image(2.5);\n",
    "    coefficient_1 = ee.Image(6);\n",
    "    coefficient_2 = ee.Image(7.5);\n",
    "    l = ee.Image(1);\n",
    "    nir = image.select(\"B5\");\n",
    "    red = image.select(\"B4\");\n",
    "    blue = image.select(\"B2\");\n",
    "    evi = image.expression(\n",
    "        \"Gain_Factor*((NIR-RED)/(NIR+C1*RED-C2*BLUE+L))\",\n",
    "        {\n",
    "            \"Gain_Factor\":gain_factor,\n",
    "            \"NIR\":nir,\n",
    "            \"RED\":red,\n",
    "            \"C1\":coefficient_1,\n",
    "            \"C2\":coefficient_2,\n",
    "            \"BLUE\":blue,\n",
    "            \"L\":l\n",
    "        }\n",
    "    )\n",
    "    return evi\n",
    "\n",
    "#Atmospherically Resistant Vegetation Index\n",
    "def ARVI(image):\n",
    "    red = image.select(\"B4\")\n",
    "    blue = image.select(\"B2\")\n",
    "    nir = image.select(\"B5\")\n",
    "    red_square = red.multiply(red)\n",
    "    arvi = image.expression(\n",
    "        \"NIR - (REDsq - BLUE)/(NIR+(REDsq-BLUE))\",{\n",
    "            \"NIR\": nir,\n",
    "            \"REDsq\": red_square,\n",
    "            \"BLUE\": blue\n",
    "        }\n",
    "    )\n",
    "    return arvi\n",
    "\n",
    "#Leaf Area Index\n",
    "def LAI(image):\n",
    "    nir = image.select(\"B5\")\n",
    "    red = image.select(\"B4\")\n",
    "    coeff1 = ee.Image(0.0305);\n",
    "    coeff2 = ee.Image(1.2640);\n",
    "    lai = image.expression(\n",
    "        \"(((NIR/RED)*COEFF1)+COEFF2)\",\n",
    "        {\n",
    "            \"NIR\":nir,\n",
    "            \"RED\":red,\n",
    "            \"COEFF1\":coeff1,\n",
    "            \"COEFF2\":coeff2\n",
    "        }\n",
    "    )\n",
    "    return lai\n",
    "\n",
    "def tasseled_cap_transformation(image):\n",
    "    #Tasseled Cap Transformation for Landsat 8 based on the \n",
    "    #scientfic work \"Derivation of a tasselled cap transformation based on Landsat 8 at-satellite reflectance\" \n",
    "    #by M.Baigab, L.Zhang, T.Shuai & Q.Tong (2014). The bands of the output image are the brightness index, \n",
    "    #greenness index and wetness index.\n",
    "    b = image.select(\"B2\", \"B3\", \"B4\", \"B5\", \"B6\", \"B7\");\n",
    "    #Coefficients are only for Landsat 8 TOA\n",
    "    brightness_coefficents= ee.Image([0.3029, 0.2786, 0.4733, 0.5599, 0.508, 0.1872])\n",
    "    greenness_coefficents= ee.Image([-0.2941, -0.243, -0.5424, 0.7276, 0.0713, -0.1608]);\n",
    "    wetness_coefficents= ee.Image([0.1511, 0.1973, 0.3283, 0.3407, -0.7117, -0.4559]);\n",
    "    fourth_coefficents= ee.Image([-0.8239, 0.0849, 0.4396, -0.058, 0.2013, -0.2773]);\n",
    "    fifth_coefficents= ee.Image([-0.3294, 0.0557, 0.1056, 0.1855, -0.4349, 0.8085]);\n",
    "    sixth_coefficents= ee.Image([0.1079, -0.9023, 0.4119, 0.0575, -0.0259, 0.0252]);\n",
    "    \n",
    "    #Calculate tasseled cap transformation\n",
    "    brightness = image.expression(\n",
    "        '(B * BRIGHTNESS)',\n",
    "        {\n",
    "            'B':b,\n",
    "            'BRIGHTNESS': brightness_coefficents\n",
    "        })\n",
    "    greenness = image.expression(\n",
    "        '(B * GREENNESS)',\n",
    "        {\n",
    "            'B':b,\n",
    "            'GREENNESS': greenness_coefficents\n",
    "        })\n",
    "    wetness = image.expression(\n",
    "        '(B * WETNESS)',\n",
    "        {\n",
    "            'B':b,\n",
    "            'WETNESS': wetness_coefficents\n",
    "        })\n",
    "    fourth = image.expression(\n",
    "        '(B * FOURTH)',\n",
    "        {\n",
    "            'B':b,\n",
    "            'FOURTH': fourth_coefficents\n",
    "        })\n",
    "    fifth = image.expression(\n",
    "        '(B * FIFTH)',\n",
    "        {\n",
    "            'B':b,\n",
    "            'FIFTH': fifth_coefficents\n",
    "        })\n",
    "    sixth = image.expression(\n",
    "        '(B * SIXTH)',\n",
    "        {\n",
    "            'B':b,\n",
    "            'SIXTH': sixth_coefficents\n",
    "        })\n",
    "    bright = brightness.reduce(ee.call(\"Reducer.sum\"));\n",
    "    green = greenness.reduce(ee.call(\"Reducer.sum\"));\n",
    "    wet = wetness.reduce(ee.call(\"Reducer.sum\"));\n",
    "    four = fourth.reduce(ee.call(\"Reducer.sum\"));\n",
    "    five = fifth.reduce(ee.call(\"Reducer.sum\"));\n",
    "    six = sixth.reduce(ee.call(\"Reducer.sum\"));\n",
    "    tasseled_cap = ee.Image(bright).addBands(green).addBands(wet).addBands(four).addBands(five).addBands(six)\n",
    "    return tasseled_cap.rename('brightness','greenness','wetness','fourth','fifth','sixth')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Derive and visualize Tasseled Cap Transformation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "tct = tasseled_cap_transformation(landsat8_collection)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"https://earthengine.googleapis.com/api/thumb?thumbid=d613d8639b93fea2a6a1acdaec28433f&token=d437a913667c9e68ddab34c3c1894b44\"/>"
      ],
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "image = tct.clip(area_of_interest)\n",
    "theGeom = area_of_interest.getInfo()['coordinates']\n",
    "thumbnail = image.getThumbUrl({'min':-1,'max':2,'size':'800','bands':['brightness','greenness','wetness'],\n",
    "                               'region':theGeom})\n",
    "Image(url=thumbnail)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Derive indices, spectral angles. Build and visualize image stack"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "ndvi = NDVI(landsat8_collection)\n",
    "sam = SAM(landsat8_collection)\n",
    "evi = EVI(landsat8_collection)\n",
    "arvi = ARVI(landsat8_collection)\n",
    "lai = LAI(landsat8_collection)\n",
    "spectral_indices_stack = ee.Image(ndvi).addBands(lai).addBands(sam).addBands(arvi).addBands(evi).addBands(tct).addBands(landsat8_collection)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"https://earthengine.googleapis.com/api/thumb?thumbid=15a39aea7c829d608ea807c7eaac0f3f&token=41d9cb3614ba5d875ab968f846ef42f0\"/>"
      ],
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "image = ndvi.clip(area_of_interest)\n",
    "theGeom = area_of_interest.getInfo()['coordinates']\n",
    "thumbnail = image.getThumbUrl({'min':-1,'max':1,'size':'800',\n",
    "                               'region':theGeom})\n",
    "Image(url=thumbnail)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define classification function\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def classification(raster_input, training_dataset,number_of_training_points, region, classification_algorithm):\n",
    "    bands = raster_input.bandNames()\n",
    "    points = ee.FeatureCollection.randomPoints(region, number_of_training_points, number_of_training_points, 1)\n",
    "    training = training_dataset.addBands(raster_input).reduceToVectors(\n",
    "        reducer='mean',\n",
    "        geometry=points,\n",
    "        geometryType='centroid',\n",
    "        scale=30,\n",
    "        crs='EPSG:4326'\n",
    "    )\n",
    "    classifier = ee.Classifier.randomForest().train(\n",
    "        features=training,\n",
    "        classProperty='label',\n",
    "        inputProperties=raster_input.bandNames(),\n",
    "    )\n",
    "    out = raster_input.classify(classifier)\n",
    "    return out"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Derive classification function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "output = classification(spectral_indices_stack, mexico_landcover_2010_landsat, 10000, area_of_interest, 'Cart')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "palette = ['5d9cd4','007e00','003c00','aaaa00','aa8000','8baa00','ffb265','00d900','aa007f','ff55ff','ff557f','ff007f','ff55ff','aaffff','00ffff','55aaff','e29700','bd7e00','966400','a2ecb1','c46200','aa5500','6d3600','00aa7f','008a65','005941','e9e9af','faff98',\n",
    "'00007f','c7c8bc','4d1009','000000','fef7ff','6daa50','3a7500','0b5923','ffaaff','ffd1fa']\n",
    "palette = ','.join(palette)\n",
    "\n",
    "# make a visualizing variable\n",
    "vis = {'min': 0, 'max': len(palette), 'palette': palette, 'region':theGeom}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Display training data of classification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"https://earthengine.googleapis.com/api/thumb?thumbid=789dbd637c1255d30e116abc65bd3669&token=faf2c48b999f387709e7e02caf9648ca\"/>"
      ],
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "image = mexico_landcover_2010_landsat.clip(area_of_interest)\n",
    "theGeom = area_of_interest.getInfo()['coordinates']\n",
    "thumbnail = image.getThumbUrl(vis)\n",
    "Image(url=thumbnail)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Display classification output\n",
    "\n",
    "Please be patient. It may take a few moments. You might have to run this cell several times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"https://earthengine.googleapis.com/api/thumb?thumbid=c184b64f683ffb1604a5d44302ca9300&token=5189c2b6909b24dfa1508f787c62ed9e\"/>"
      ],
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "image = output.clip(area_of_interest)\n",
    "theGeom = area_of_interest.getInfo()['coordinates']\n",
    "thumbnail = image.getThumbUrl(vis)\n",
    "Image(url=thumbnail)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
